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Abstract. The interactions between climate and the environment are highly 
complex. Due to this complexity, process-based models are often preferred 
to estimate the net magnitude and directionality of interactions in the Earth 
System. However, these models are based on simplihcations of our under¬ 
standing of nature, thus are unavoidably imperfect. Conversely, observation- 
based data of climatic and environmental variables are becoming increasingly 
accessible over large scales due to the progress of space-borne sensing tech¬ 
nologies and data-assimilation techniques. Albeit uncertain, these data en¬ 
able the possibility to start unraveling complex multivariable, multiscale re¬ 
lationships if the appropriate statistical methods are applied. 

Here, we investigate the potential of the wavelet cross-correlation method 
as a tool for identifying multiscale interactions, feedback and regime shifts 
in geophysical systems. The ability of wavelet cross-correlation to resolve the 
fast and slow components of coupled systems is tested on synthetic data of 
known directionality, and then applied to observations to study one of the 
most critical interactions between land and atmosphere: the coupling between 
soil moisture and near-ground air temperature. Results show that our method 
is not only able to capture the dynamics of the soil moisture-temperature 
coupling over a wide range of temporal scales (from days to several months) 
and climatic regimes (from wet to dry), but also to consistently identify the 
magnitude and directionality of the coupling. Consequently, wavelet cross¬ 
correlations are presented as a promising tool for the study of multiscale in- 
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teractions, with the potential of being extended to the analysis of causal re¬ 
lationships in the Earth system. 
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1. Introduction 

Natural systems are often characterized by complex interactions, mainly originating in 
the overlap of dynamical processes acting at very diverse temporal and spatial scales. Ex¬ 
amples of this multiscale dynamics can be found in several branches of geophysics. These 
include climate and the hydrological cycle, whose different components interact and syn¬ 
chronize over a wide range of scales and patterns [Lovejoy and Schertzer, 2013; Tsonis 
and Eisner, 2007], and ecological systems, where resilience and evolntion are mainly de¬ 
termined by cooperation and connectivity [Sole et ai, 1998; Moilanen and Nieminen, 
2002; Cowen et ai, 2006]. In a similar framework it is logical to expect that also the 
conplings among the different components of the system - and between the system and 
the surronnding - take place at mnltiple scales and across them. Mnltiscale interactions 
have recently received extensive attention in the literatnre, and have been proposed as a 
mechanism for the triggering of extreme events [Miralles et ai, 2014; Peters, Debra P. C. 
et ai, 2004; Rajfa et ai, 2008], abrupt regime transitions [Okin et ai, 2009; Peters et ai, 
2007] and patterns formation [Scanlon et ai, 2007; Guttal and Jayaprakash, 2009]. 

Examples of this increasing interest for mnltiscale and cross-scale interactions can be 
fonnd in ecology [Allen and Plolling, 2013; Moritz et ai, 2005; Cash et ai, 2006; Peters 
et ai, 2007; Raffa et ai, 2008; Scanlon et ai, 2007; Thrush et ai, 2013; ITerner et ai, 2014] 
and climate dynamics [Plolbrook et ai, 2014; Debra P. C. Peters et ai, 2007; Molini et ai, 
2010a; Okin et ai, 2009; Rial et ai, 2004], but also in helds other than geosciences such 
as network morphology [Odor, 2013; Pastor-Satorras and Vespignani, 2001] and econo¬ 
metrics [Nikkinen et ai, 2011]. Most of these stndies are based on minimalist models of 


DRAFT 


February 19, 2015, l:27ain 


DRAFT 



X - 6 


CASAGRANDE ET AL.: MULTISCALE COUPLING IN THE EARTH SYSTEM 


interaction across scales [Allen and Holling, 2002; Peters, Debra P. C. et ai, 2004; Peters 
et ai, 2007], or - when some kind of data-driven approach is attempted - on classical 
scaling statistics, more able to resolve the scale-dependent structure of the considered 
processes, rather than coupling strength and directionality across scales. 

Following this approach, the scaling properties of a wide number of dynamical processes 
including turbulent flows [Frisch and Kolmogorov, 1995], atmospheric tracers such as pre¬ 
cipitation [Molini et ai, 2010b; Veneziano and Lepore, 2012; Lovejoy and Schertzer, 2013], 
ecosystems patterns and organization [Wu, 2006; Nagelkerken, 2009], social networks [Szell 
et ai, 2010], urban growth and development [hammer et ai, 2006; Chen and Zhou, 2008; 
Bettencourt et ai, 2007; Pumain, 2004], and economic systems [Lux and Marchesi, 1999; 
Mandelbrot and Stewart, 1998; Mandelbrot, 1997] have been intensively investigated in 
the last decades through a variety of scaling metrics. In contrast, considerably less atten¬ 
tion has been devoted to the analysis of couplings and feedbacks taking place at different 
scales and across them [Dhamala et ai, 2008a; Molini et ai, 2010a; Schmitt and Chainais, 
2007], to the development of ad-hoc statistics to analyze the temporal evolution of such 
couplings [Li and Nozaki, 1997; Ngae et ai, 1998; Mizuno-Matsumoto et ai, 2001; Sello 
and Bellazzini, 2000; Salvetti et ai, 1999; Onorato et ai, 1997; Lungarella et ah, 2007; 
Shirazi et ai, 2013], and to the investigation of the spectral features of systems displaying 
strong connectivity across inter-annual, seasonal and sub-seasonal scales [Torrence and 
Compo, 1998; Torrence and Webster, 1999]. 

In this study we explore the potential of a simple local-coupling metric, the wavelet 
cross-correlation, for identifying and assessing linear and intermittent interactions across 
different temporal scales in geophysical systems. Time-domain correlations are routinely 
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applied to multivariate geophysical time-series to identify phase consistency among the 
variables, but they fails when the variability is dominated by periodicity or when the phase 
consistency is intermittent. Harmonics-based statistics eqnivalent to the time-domain 
correlation, like spectral coherence and cross-spectra can partially £11 this gap, being able 
to separate the phase consistency over different persistent freqnencies. However, also these 
spectral statistics are limited in captnring interactions that may be transient in time or 
may display different directionality across scales. 

Conversely, wavelet cross-correlation can be inferred throngh the wavelet decomposition 
of observed signals, and widens the classical concept of multiscale information flow within 
random mnltiplicative processes [Arneodo et ai, 1998]. What is particularly appealing 
in this measnre is its ability to decompose linear correlations in scale and time, preserv¬ 
ing simultaneously the total correlation of the system [Daubechies, 1992]. Based on the 
whitening properties of wavelet filters [Mallat, 2008; Percival, 1999], the wavelet coeffi¬ 
cients describing the bi-variate process {Xt, Yt} at each temporal scale s can be treated as 
realizations of the jointly Ganssian random process {Wx(t, s), Wy (t, s)} and the sample 
pairs {wx(ti, s),WY(ti, s)) and {wx(tj, s),WY(tj, s)) can be considered mntnally indepen¬ 
dent for i ^ j. However, if the joint-Gaussianity assnmption can in practice be relaxed 
when we test nnll cross-correlations, the null auto-correlation hypothesis snffers from a 
nnmber of limitations - mostly depending on the degree of wavelet overlapping at the 
analyzed scale - which need to be considered. 

The essential theoretical backgronnd on wavelet cross-correlation and its significance 
testing are discussed in section 2, together with some caveats abont the inference of di¬ 
rectional conplings from spectral statistics. Section 3 is devoted to test the performance 
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of the wavelet cross-correlation in assessing scale-by-scale interaction in synthetic auto¬ 
regressive systems of known directional coupling and increasing complexity. Following the 
test on synthetic data, section 4 deals with the multiscale nature of land-atmosphere inter¬ 
actions through the example of the air temperature-soil moisture coupling and feedback 
[Seneviratne et ai, 2010, 2006; Miralles et ai, 2012; Orlowsky and Seneviratne, 2010]. Fi¬ 
nally, a further discussion of wavelet cross-correlation strength and criticality is provided 
in section 5, together with concluding remarks and future developments. 


2. Methods and Background 

Wavelet cross-correlations are obtained through the combined use of wavelet hltering 
[Mallat, 2008; Daubechies, 1992] and classic linear coupling measures. 

In the time domain, under stationary and ergodic assumptions, the simplest and most 
adopted measure of linear coupling between the trajectories of a bi-variate real-valued 
stochastic process {Xt,Yt} is the cross-correlation function pxvij) dehned as: 

7au(u) 


Pxy{t) = 




( 1 ) 


with r = 0,..., A; representing the time lag (i.e. temporal asymmetry) between the process 
trajectories, 7 xu(t) = E[XtYtJ^r\ their non-centered covariance, and ay the variances 
of {W|'} and {'lYt} respectively. Assuming that causes precede effects in time, it is com¬ 
mon practice to associate the presence of signihcant correlations at non-zero lags with 
causal asymmetric coupling. However, when the analyzed signals display multi-scaling, 
non-stationarity and periodicity - or simply some form of oscillatory behavior at differ¬ 
ent frequencies - causality can hardly be inferred from classic lagged cross-correlations 
[see Granger, 1969, for a detailed discussion of the relationship between causality, co- 
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integration and correlation]. In fact, it is important to note that even if the analyzed 
signals are expected to display different coupling strengths and synchronization patterns 
over a wide range of temporal scales - like often the case in climate and other geophys¬ 
ical systems - what we dually “see” through the estimation of pxvit-it + u) is only the 
aggregated effect of these multiscale interactions, that does not necessary reproduce the 
actual directionality of the coupling at a specidc scale s. 

As an example, interactions between land and atmosphere such as the soil moisture-air 
temperature coupling discussed in section 4, are characterized by strong seasonal syn¬ 
chronization effects that can partially mask the directionality of due-scale interactions. 
However, high-frequency components - through vertical duxes of energy and water acting 
at extremely localized spatio-temporal scales - still play an important role in trigger¬ 
ing and sustaining those interactions. Decomposing the correlation into its scale-by-scale 
components can hence shade some light on the diderent dynamical processes character¬ 
izing the observed coupling. This can be achieved by simply considering the analogue 
of lagged correlation in the wavelet domain [Li and Nozaki, 1997; Turbelin et ai, 2008], 
i.e. wavelet eross-correlation. Unlike integral statistics such as wavelet co-spectra and 
coherence [Grinsted et ai, 2004; Maraun and Kurths, 2004; Torrenee and Compo, 1998], 
the wavelet cross-correlation is based on the “a priori” decomposition of the bi-variate sig¬ 
nals through wavelet band-pass dltering and on a direct inference of scale-by-scale linear 
correlations from the resulting time series of sample coefdcients WXi{t,s). 

Wavelet cross-correlations were drst adopted in a number of experimental studies on 
multiscale interactions in turbulent dows and mixing [Li and Nozaki, 1997; Ngae et ai, 
1998; Sello and Bellazzini, 2000; Salvetti et ai, 1999; Onorato et ai, 1997], mainly based 
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on the continuous wavelet transform (CWT). Following this bulk of work, Whitcher and 
coauthors [Whitcher and Jensen, 2000; Whitcher et ai, 2000] proposed an expression for 
wavelet cross-correlations and their conhdence interval based on the maximum overlap 
discrete wavelet transform (MODWT), a variant without sub-sampling of the orthonor¬ 
mal discrete wavelet transform (DWT) [Foufoula-Georgiou and Kumar, 1994; Percival 
and Mofjeld, 1997]. The use of MODWT instead of DWT was dictated by the lack of 
translation invariance in the latter, strongly impacting the hnal lag-resolution of wavelet 
cross-correlations. Through the MODWT, it is also possible to reduce the effects of re¬ 
dundancy, partially preserving invariance in the translation [Gengay et ai, 2001]. 

However, the trade-off between lag-resolution and redundancy should be carefully eval¬ 
uated depending on the specihc application. In this study, based on the wide range of 
temporal scales over which geophysical processes evolve, and on the necessity of retain¬ 
ing as much information as possible about asymmetry in coupling, we opted for using 
the CWT together with complex analyzing wavelets, known to preserve these properties 
[Mallat, 2008]. Also redundancy of CWT cannot be really considered a disadvantage here, 
until we are not concerned with compression and the effects of auto-correlation at large 
scales. Rather it can often become an advantage since the redundancy of CWT allows for 
a better visualization of correlation patterns across scales [Growley, 2007]. 

2.1. Continuous Wavelet Filtering: Generalities 

Let x(t) denote a sample trajectory of the finite energy random process {Xt}. Then 
via the CWT, we can decompose x{t) into a set of hnite basis functions, representing 
its variability at different scales and instants in time [Mallat, 2008; Daubechies, 1992; 
Torrence and Gompo, 1998]. 
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The coefficients wx{u,s) of the CWT are obtained by decomposing x(t) over dilated 

and translated wavelet functions ip G L^(M) of zero average and ||'^|| = 1: 

/ + 00 

x{t)^ 

VS 

where s is the wavelet scale (inverse of the pseudo-frequency), u is the translation along 
the time axis, ip* (■) indicates the complex conjugate of the wavelet basis function and {•, •) 
is the inner product. Since each trajectory of {Xt} is a sample function from a collection of 
random variables the wavelet hltering process is asymptotically equivalent to 

decomposing {Xt} in a hnite number of stochastic processes {VVx(m, <s)} whose realizations 
are the wx{u, s). In addition, it is easy to prove that the wavelet transform is equivalent 
to a convolution with dilated band-pass hlters [Mallat, 2008, p. 79]. 

Thus, we can rewrite equation (2) as wx{u, s) = x-kipsiu), with ipsiu) = 1/^/s 1 p* {—t/s) 
and its Fourier transform given by ipg{u}) = \Psip* {su). It follows that since ip is zero- 
average, V(0) is also zero and ip becomes the transfer function of a band-pass hlter. The 
CWT extends the benehts of Fourier analysis to observations involving transients and 
non-stationarities, allowing the estimation of local spectral density metrics such as the 
wavelet scalogram: 

Wx{t, s) = \wx{t, s)\^ , (3) 

and cross-scalogram; 

hFxv(t, s) = w*x{t, s)wY{t, s), (4) 



routinely used as a measure of coupling across scales. However, these and other similar 
metrics such as the wavelet co-spectrum and the wavelet coherence, are not able to provide 
any information about the temporal asymmetry of couplings at different scales, unlike 
wavelet cross-correlation metrics. 
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2.2. Multiscale Interactions and Wavelet Cross-correlations 

We can now define the wavelet covariance of the random process {Xf, 1^} at lag r and 
scale s, as: 

7xy{s,t) = E[Wx{t,s)WY{t + T,s)] (5) 


This can be alternatively expressed as [Li and Nozaki, 1997]: 


'*+00 


7xu(s,'r) = — / SxY{uj)\ij{suj) 




where Sxy{^) is the co-spectrum of the two signals and 


rxY+,s) = Sxy{uj) 'Ipisu) 


( 6 ) 


( 7 ) 


is the local wavelet co-spectrum function. Therefore 7 xY("S,r) is the inverse Fourier trans¬ 
form of the local co-spectrum, that integrated across scales gives the classical covariance 
among the signals 7 xy(t)- Although, if the analyzing wavelet tjj is complex also 7x^(5, r) 
is a complex function and can be decomposed into a real part 3? ( 7 x^( 5 , r)) and an imag¬ 
inary part A ( 7 xy(<s, r)), bearing information about the strength and the phase of the 
correlations. It is important to note that in this last case the conservation of 7 xy(t) 
across scales holds only for the real part of the coefficients [Daubechies, 1992], so that the 
corresponding wavelet cross-correlation pxy{s,t) is most commonly estimated as [Li and 
Nozaki, 1997]: 


Pxy{s,t) = 


3? (7xY(-s,r)) 


( 8 ) 


v/3f?(a+(s))3fJ(a+(s))’ 
where the represent the variance of the i-th variable coefficients at scale s. This form 

of wavelet cross-correlation ranges between —1 and 1 like its time domain counterpart, and 
represents a simple local decomposition of the cross-correlation function in time. Since the 
CO- variance of a bi-variate complex-valued random process {Wx(t, s), Wy (t, s)} is given 
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by: 


7 = (Wf + WlWp) + i - W^W§) , 


(9) 


it is easy to show through symmetry considerations that 


S ( 7 xy (s, t) ) = 2B [S ( Wx (« , s)) S ( Wy (« + T, s))]. 


( 10 ) 


Therefore, if an analytic wavelet is used in the decomposition, the asymmetry in the 
coupling at different scales mainly derives from temporal shifts among the real parts of 
the coefficients, rather than from a direct wavelet phase estimation. Figure 1 provides 
a conceptual representation of wavelet correlation patterns corresponding to X driving 
Y (d), instantaneous coupling (e) and Y driving X (f) and the associated geometry of 
sample coefficients at generic scale sq (a-c). Here the forcing direction is assumed to be 
homogeneous across scales like in a red-noise dominated process lacking of characteristic 
forcing scales. 

There are a number of alternative formulations for pxy{s,t), such as the one recently 
proposed by Shirazi et al. [2013], and based on the amplitudes instead of the real part of 
the coefficients. This neglects the information about temporal asymmetries to focus only 
on the amplitude of correlations across diverse scales. Sello and Bellazzini [2000] also 
proposed to estimate scale-by-scale couplings in the form of a local wavelet coherence: 

2||7Ay(r,.)f 

+ lkV(s)||^’ 

This metric however, varies between 0 and 1 due to the absolute operator at the nomi¬ 
nator, thus it cannot provide information about the sign of the coupling. 

In the following, pxy{s,t) is estimated based on the expression in equation (8), which 
is also the one most directly connected with linear couplings in the time domain. The 
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corresponding estimator of Pxy{s,t) is then given by: 


^ (E u!x (b s) w'y (i + r, g)) 

^ (E Kx (b s)]^) ^ (E ku (b s)]^) 


(12) 


with w' {i, s) = w {i, s) — Wg and Wg being the long-term average of coefficients at scale s. 


2.3. Complex Wavelet Kernels 

Different analyzing wavelets can lead to diverse localization effects in frequency and/or 
in time, that can become crucial for the identihcation of scale-dependent directional in¬ 
teractions. For this reason, in the following we estimate sample correlations rxy(s,r) by 
using two analytic wavelets with different localization in time and scale, i.e. the Morlet 
and the Paul wavelets. The Morlet wavelet is dehned as: 


^morl (h) = TT 




(13) 


while the Paul wavelet is given by: 

om-m.™ I 

Vu.O)= . ,, ' n - (14) 

v'ir(2m)! 

where r] is the non-dimensional time parameter, uq the central frequency of the Morlet 
wavelet and m the order parameter of the Paul wavelet. 

The Morlet wavelets are more localized (displaying higher resolution) in frequency than 
in time, while the Paul wavelets display a higher temporal localization [Torrence and 
Compo, 1998]. This comparison allows us to address the effects of time/scale resolution 
on the estimation of the wavelet cross correlations. Following the classic work of Torrence 
and Compo [1998] we use uq = 6 and m = 4 to obtain an effective trade-off in the 
time and frequency resolution of the decomposition. When working with the Morlet 
wavelet, the choice of a central frequency Wq = 6 is further justihed by the fact that 


Morlet is not a proper wavelet, as its integral is not zero. However, for uq larger than 
DRAFT February 19, 2015, l:27ain DRAFT 



CASAGRANDE ET AL.: MULTISCALE COUPLING IN THE EARTH SYSTEM 


X - 15 


5, the integral becomes small enough to ensure the numerical applicability of the Morlet 
kernel [Kumar and Foufoula-Georgiou, 1997]. As complex functions, both the Morlet 
and Paul wavelets are suitable to study oscillatory time series, and have been extensively 
used in geosciences [Foufoula-Georgiou and Kumar, 1994; Torrence and Gompo, 1998]. In 
addition, both the Morlet and Paul wavelets are symmetric, allowing for a non-distorted 
estimation of temporal shifts. This is an important property for robust and reliable 
analysis of directional couplings. At the same time, we have to keep in mind that being 
non-orthogonal, their CWT can be affected by an overlapping of sub-frequency bands and 
consequent redundancies in the decomposition of the analyzed signal. 

2.4. Significance Test 

Similar to their time-domain counterparts, multiscale correlations can be tested for sig- 
nihcance based on a scale-dependent approach. In the following, we introduce a simple 
signihcance test for wavelet cross-correlations based on the assumption that the coef- 
hcients {wxiti, s),WY(ti, s)) are sample pairs drawn from the jointly Gaussian random 
process {Wx(t, s), Wy(t, s)}, and that rxv(s,r) hereafter for simplicity) is the cor¬ 
responding sample statistic for the correlation strength and direction at scale s. Such 
an assumption originates from the whitening properties of the wavelet transform [Mallat, 
2008; Percival, 1999], and provides us with the necessary machinery to test ^ against the 
null correlation hypothesis and construct approximate statistical conhdence intervals. We 
show that in case of null correlation the joint-normality condition can be relaxed owing 
to the asymptotic properties of the sample correlation distribution in ps^r = 0 [Johnson 
et ai, 1994]. 
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In general, the probability distribution of the sample correlation fRlvs^r) is cumbersome 
to derive in a closed-form. One of the simplest expressions is due to Fisher [Fisher, 1915] 
and is given by: 


/R(^s,r) 


n —1 

-Pjr) ^ 
ttF (n — 2 ) 


— ) 

S,T J 


n —4 
2 


^ I COS ^{-rs,rPs,T) I 

d{rs,rps,r) I y^(l J J 

where r(-) is the Gamma function and n the sample size. Equation (15), although written 
in terms of elementary functions, is still too complex to be explicitly used in testing. 
However, for ps^r = 0 it reduces to the null pdf proposed in 1908 by Student [Kendall and 
Stuart, 1945; Johnson et ah, 1994]: 


fnirs,-. 




n —4 
2 


r(¥) 

r(i) 




n —4 
2 


(16) 


where B{ ) is the Beta function. Contrary to sample distributions for p 7 ^ 0 - known to 
be markedly asymmetric ~ the probability density function in equation (16) is symmetric 
around 0 and its derived distribution for 


t T 


V(n-2), 


(1 -^s,r) 

reduces to a t-Student with (n — 2 ) degrees of freedom: 


(17) 


frits, t) — 


1 + 


n —1 

2 


(18) 


n-2. 

Therefore, once dehned a significance level a, equations (17) and (18) can be used to 
test wavelet cross-correlations rs,T against the hypothesis Hq : ps,r = 0 , whose two-tailed 
rejection region lies outside [—ta/2, n-2, ta/2, n-2] • An alternative approach, often adopted 


DRAFT 


February 19, 2015, l:27ain 


DRAFT 



CASAGRANDE ET AL.: MULTISCALE COUPLING IN THE EARTH SYSTEM 


X - 17 


in the literature to test the hypothesis p > p' for p' ^ 0, relies on a variance-equalizing 
transformation known as Fisher ; 2 -transformation; 

Zs^r = tanh" ^ log Q ^ ^ , (19) 

where Z is approximately normally distributed with mean = (1/2) ln[(l-|-ps,T)/(l—P s,t)] 
and variance a1 = l/(n — 3). However in this study we limit the approach to test the 
absence of correlation across scales. 

This choice is motivated by the superior robustness of the estimator in p = 0, and 
to the possibility of widening the usage of null distribution in equation (16) to non jointly- 
Gaussian samples [Johnson et ai, 1994], In fact, while the expression in equation (16) 
can be seen as an exact representation of fn in p = 0, Fisher’s transformation only 
represents an approximation of fn away from 0. Consequently, assuming that linearity 
and homoscedasticity conditions hold, the //? in p = 0 always coincides with the null 
distribution of a jointly-normal process, while the non null distribution of r is robust only 
under additional kurtosis constrains [see Johnson et ai, 1994, p.582]. It is also worth 
noting that the transformation in equation (18) is generally assumed to be less sensitive 
to violations of the normality assumption [Edgell and Noon, 1984; Havlicek and Peterson, 
1977], i.e. deviations of {Wx(f, s), Wy (t, s)} from a jointly Gaussian distribution should 
not impact the test for null correlation across scales. 

In geosciences it is also a common practice to test spectral statistics - such as wavelet 
co-spectra and coherence - against alternative null hypotheses based on the background 
noise of the underlying process [Torrence and Webster, 1999]. The statistical signihcance 
at the desired signihcance level a is therefore obtained numerically via Monte Garlo simu¬ 
lations. For example, in the case of a red background noise, the null hypothesis is obtained 
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by Monte Carlo replicates [Torrence and Webster, 1999] through the following steps: (a) a 
hrst order autoregressive model of the background noise is htted to the observed data; (b) 
a set of surrogates is generated from the htted model and (c) a suitable conhdence interval 
is computed producing the desired quantile bounds. However, this numerical approach 
has been shown to lead to spurious signihcant correlations in case of reduced samples [Ma- 
raun and Kurths, 2004; Maraun et ai, 2007]. Also, it implies strong assumptions on the 
background noise of the observed processes. 

3. Wavelet Cross-correlation from Auto-regressive Systems with Pseudo- 
periodic Features 

In this section we test the ability of the wavelet cross-correlation rg^r to capture multi¬ 
scale interactions in synthetic processes of known directional coupling. Large ensembles 
of numerically generated time series are used to understand whether rs,r can be used as 
an efficient estimator of scale-by-scale coupling strength and directionality for systems of 
increasing complexity. 

3.1. Coupling in a First-Order Vector Auto-regressive System 

The hrst case study we consider is a simple hrst order vector auto-regressive model 
(VAR(l)) in the form: 

x{t) _ 0.80 Cl x{t — 1) e(t) 
y{t)_ ~ _C -2 0.70 

where Ci and C 2 are coupling parameters dehning the strength and directionality of the 
interaction between x and y, and e{t) and ^{t) are uncorrelated noise terms with zero mean 
and unitary variance. Figure 2a shows a realization of the VAR(l) model in equation (20), 
where the x and y sub-spaces were linked in a unidirectional way by imposing a coupling 
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coefficient Ci = —0.4 {y ^ x), and a null feedback from x to y (C 2 = 0). Panel b of the 
same figure reports the corresponding theoretical power spectra for the x (blue dashed 
line) and y (green solid line) autoregressive sub-spaces, whose Lorentzian decay is typified 
by the absence of characteristic scales of autocorrelation. VAR(l) models, even though 
simple in construct, still represent a valuable benchmark for Ts^t- Red noise, or more in 
general l//“ noise spectral decay is in fact a characteristic feature of many geophysical 
systems [Agnew, 1992; Keshner, 1982; Muzzy et ai, 2011]. In these models the strength 
and directionality of the coupling can be easily - and intuitively ~ tuned. Also, the absence 
of a characteristic scale for the coupling is reflected in “cascade-like” r-g^T patterns similar 
to the ones sketched in Figure 1. 

Figure 3a-d shows the computed from both single and ensemble-realizations of the 
VAR(l) model in equation (20) for different values of the coupling parameters Ci and C 2 . 
Each ensemble includes 100 realizations of sample size n = 4096, independently generated 
starting from random initial conditions. The wavelet cross-correlation is estimated as 
an ensemble average over all the 100 realizations based on a Paul mother wavelet of order 
m = 4, and thus represented as a function of temporal asymmetry r on the abscissas and 
scale s on the ordinates. Figure 3a depicts the scale-by-scale correlation patterns resulting 
from imposing a strong negative coupling from x to y {C 2 = —0.9) and null feedback from 
y to X {Cl = 0), while Figure 3b shows a case with opposite and weaker coupling, i.e. 
Cl = —0.4 and C 2 = 0. The considered VAR(l) systems are therefore both negatively 
coupled but with different strengths, and feedback effects are not considered at this stage. 
Correlations below the a = 99% signihcance level are masked in white. 
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It is interesting to note how in a vector system such as the one in equation (20), char¬ 
acterized by red-noise-like spectral features and “short-memory”, the coupling propa¬ 
gates along all the scales in a homogeneous way. To better highlight possible asym¬ 
metries in the coupling at different temporal scales s and the variability in asymmetry 
across the diverse simulations of the ensemble, we also extracted the minimum correlation 
f'min = Hlin {rs,T}^^_f. for each single realization and derived the ensemble minimum mean 
fmin = and standard deviation The focus is here on the minimum of 

since the two synthetic systems are negatively coupled. However similar considerations 
on the scale-by-scale maximum of could be done in the case of a positively correlated 
process. As discussed next, fmin and provide a rough estimate of the asymmetry and 
variability of peak correlations in the analyzed system. 

In Figure 3a-b,d, fmin (black empty circles) and its conhdence interval at the 99% 
signihcance level (bounded by the blue and red solid lines) were overlapped to the 
diagram to show how the peak correlation moves from one side to the other of the 0 
bisecting line when the coupling direction is inverted. Top panels in Figure 3 also show 
the fast - quasi exponential - oscillatory decay of at different scales s, typical of the 
lagged-cross correlation of an absolutely integrable signal. Depending on the scale, this 
oscillation decays below the a = 99% signihcance level more (hue-scales) or less (coarse 
scales) rapidly, consistently with the fact that the memory of the process is in general 
weaker at hner scales. We also test the ability of for identifying null-coupling across- 
scales by posing Ci = C 2 = 0 in equation (20). The wavelet cross-correlation patterns 
obtained from a single-realization of the uncoupled system are shown in Figure 3c, while 
panel d depicts the corresponding ensemble statistic. As evident the ensemble r-g^T in 
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Figure 3d never passes the test for the non-null correlation and can robustly detect the 
lack of coupling across scales. Also, f min oscillates around r = 0 with strong variability 
at all scales. 

Figure 3e hnally shows the distribution of the minimum value of the at two selected 
scales, i.e. 10 and 120 samples, in the case when both coupling and direction are set 
to vary. Also here, the conhdence intervals of f min are computed from an ensemble of 
100 realizations of the VAR(l) at the level of 99%. Panel e shows that for a sufficient 
strong value of the coupling, f min can effectively capture the switch in directionality of the 
coupling at both sample scales. At larger scale (120 samples) however, while the value of 
fmin is higher than for the smaller scale (10 samples), the variability is also higher. For 
any value of the coupling parameter falling in the interval [Ci > —0.2, C 2 > —0.2] the 
estimated range of variability of fmin ends up crossing the zero lag axis and thus does not 
provide a valid estimate of the directionality of the process. Therefore, the weaker the 
coupling, the larger the uncertainty on the actual directionality. In addition, uncertainty 
increases with scale due to the higher redundancy of wavelet coefficients. 

3.2. Coupling in a Second-Order Vector Auto-regressive System with 
Localized Pseudo-periodicity 

In this section we test the capability of Tg^r in detecting couplings that are well-localized 
in frequency - i.e. occurring at a specihc time-scale. The ability to resolve such localized 
correlations is important to the study of environmental systems, which often display strong 
sub-seasonal, seasonal and inter-annual oscillations. To do so, we consider the following 
second order vector autoregressive model (VAR(2)) where a pseudo-periodic coupling from 
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1 / to x is imposed by choosing roots close to the unit circle: 
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A similar system is used in Dhamala et al. [2008b] in the context of spectral Granger 
causality metrics. Given the localization in frequency of the simulated coupling, the 
system in equation (21) is additionally used to demonstrate the role of wavelet localization 
in the efficient estimation of Ts^t- 

Top panels in Figure 4 show as estimated from an ensemble of 100 realizations 
of equation (21), and by respectively using the Morlet (a) and Paul (b) wavelets. Both 
plots only display values above the a = 99% signihcance level, and they both show the 
conhdence interval (always 99%) for f min averaged over the 100 realizations. Panel c 
reports the theoretical power spectra of each of the x and y autoregressive sub-spaces 
displaying a clear periodicity at 1/5 of the cycle. Figure 4a-b clearly show that the r^ ,- 
obtained through the Morlet wavelet better identihes the frequency peak of the coupling 
in agreement with the characteristics of the kernel, that is by construction, more localized 
in frequency. In contrast, the frequency peak in the Paul results more spread out, 
especially in the range of lower scales, but better resolved in time. 

Nonetheless, both the Paul and the Morlet Ts^t are able to correctly detect the direc¬ 
tionality {y ^ x), as evident in the left-asymmetry of fmin. It is also worth noting that 
the variability of fmin is lower where the coupling is stronger (around the frequency peak), 
gradually increasing with the scale. This is mainly due to the size and relative higher over¬ 
lapping of the wavelets at larger temporal scales. In summary, both the Morlet and the 
Paul wavelet cross-correlations are able to extract the essential features of the coupling. 
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The choice between the two different wavelets should be made based on the application 
and the characteristics of the signals under investigation. 

3.3. Fast/Slow Dynamics Separation and Feedback 

In section 1 we argue that the evolution of geophysical systems often results from the 
interaction of diverse dynamical scales. A simple example can be provided by a coupling- 
feedback system in which the forcing is acting at shorter temporal scales than the response 
- i.e. a fast/slow dynamical system. In such a case, classic correlation analysis in the time 
domain fails unless the strength of the coupling in one of the two directions is dominant. 
The wavelet cross-correlation, in contrast, can still resolve the two components of the 
coupling, their asymmetry and characteristic scales. Let us consider a modification of 
the VAR(2) system in equation (21) in which x is driving y (negative correlation) at a 
frequency of 1/15 of cycle {fsiow) and y is forcing x (positive correlation) around the 1/6 
of cycle (// 

asi) * 

Ttl ^ [1.73 O.OOl [Tt-il r-0.90 -O.lOl [^^-21 [etl 

_yt\ ~ [O.OO 0.85J [yt.,\ + [ 0.30 -O.OSJ [yt. 2 \ ^ ^ 

The corresponding estimated by using both the Morlet and the Paul wavelets is 
shown in Figure 5, together with the normalized spectra of the two variables. It is evident 
that also in this case, the r* ^ is able to resolve the scales at which the coupling actually 
take place, as well as the opposite asymmetry of the correlation at these scales. In this 
case the f r^^r, is only representative of the coupling from x to y due to the negative 
linear correlation among the two variable at fsiow However, mainly as a consequence 
of the strong periodicity of the coefficients at both the forcing and feedback scales, the 
localization of fmin is also at some extent representative of the asymmetry in the coupling 
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from y to x. This is overall an interesting result considering that interactions characterized 
by different degrees of “memory” and fluctuation-response relaxation (FRR) effects are 
ubiquitous in geosciences [Lacorata and Vulpiani, 2007; Leith, 1975]. 

4. Multiscale Interactions in the Soil Moisture-Temperature Coupling 

Land-atmosphere interactions, their strength and directionality, are one of the main 
sources of uncertainty in current climate models with strong implications for the accurate 
assessment of future climate change impacts at regional scales [see e.g. Seneviratne et ai, 
2010]. Besides the scarcity of direct observations of the states and fluxes across the land- 
atmosphere continuum, major uncertainties originate from the inherent complexity in 
the way these variables interact, the multiscale character of these interactions, and the 
existence of critical tipping points in water and energy availability that may trigger regime 
transitions. In this last section, we apply the wavelet cross-correlation analysis to a classic 
form of interaction between land and atmosphere: the coupling between soil moisture {6) 
and near-surface air temperature (T). 

The objective is to isolate the different components of the coupling across a wide range 
of temporal scales (from fine weather scales to seasonal to inter-annual) and considering 
different lag-times between the variables. These two variables {6, T) are mainly related 
through the process of latent heat flux. They are a priori negatively correlated; however, 
their coupling can occur in both directions: (a) T may regulate 6 via the drying of the 
soil due to evaporation, or (b) 9 may regulate T due to evaporative cooling [Seneviratne 
et ai, 2010; Miralles et ai, 2012; Mueller and Seneviratne, 2012]. The latter feedback 
of 0 on T occurs in regions that are water-limited and it has been referred as a reason 
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why droughts and heatwaves coexist, or follow one another [Quesada et ai, 2012; Mimlles 
et al., 2014], 

Nonetheless, the mechanisms through which different dynamical scales contribute to 
the onset and persistence of the 6-T coupling and feedback remain unclear to a large 
extent [Orlowsky and Seneviratne, 2010]. In this study we use estimates of daily max¬ 
imum air temperature obtained from the sub daily screen level (2 m) temperatures of 
ERA-Interim (http://apps.ecmwf.int/datasets/data/interim full daily/) the most recent 
climate reanalysis product of the European Center for Medium Range Weather Forecasts 
(ECMWF) [Dee et al, 2011]. We consider temperature fields for the period 1980-2011 over 
a global grid with a spatial resolution of 0.75 degrees, corresponding to the native reduced 
Gaussian grid of ERA [Berrisford et al, 2011; Simmons et al, 2014]. Temporal resolution 
of the assimilated and predicted fields (i.e. the analysis fields) we use here is 6 hours. 
We extract and analyze maximum daily temperatures as they are most strongly impacted 
by soil moisture deficits and evaporative cooling effects [Orlowsky and Seneviratne, 2010; 
Mueller and Seneviratne, 2012]. It is important to note that ERA-Interim temperature is 
constrained by observations through a complex 4-D assimilation process [Dee et al, 2011]. 
As a consequence the quality of the analyzed fields is strongly dependent on the density of 
the station network in the region under consideration [see Simmons, 2011]. However, the 
main advantage of reanalysis over direct observations and classical interpolation schemes 
is their ability to combine observations with a physical model of the atmosphere able 
to produce physically coherent high-resolution fields and propagate information to areas 
with poor observational coverage. For 6 we use an independent data source, the global 
daily root-zone soil moisture estimates from GLEAM (Global Land-surface Evaporation: 


DRAFT 


February 19, 2015, l:27ain 


DRAFT 



X - 26 


CASAGRANDE ET AL.: MULTISCALE COUPLING IN THE EARTH SYSTEM 


the Amsterdam Methodology - http;//foofoo.ugent.be/satex/GLEAM/ as described in 
Miralles et al. [2012, 2011a, b, 2013], also for the period 1980-2011. GLEAM is a set 
of algorithms designed to retrieve information on evaporation from current satellite ob¬ 
servations of hydro-climatic variables [Miralles et al, 2013]. In GLEAM, 6 is derived 
at daily timescales through the assimilation of satellite soil moisture observations into a 
multi-layer running water balance that reproduces the inhltration of rainfall through the 
vertical soil prohle [Owe et al, 2008]. 

Figure 6 shows the results of the wavelet cross-correlation analysis for different geograph¬ 
ical locations, spanning a wide range of diverse climatic regimes. From top to bottom the 
panels are ordered by decreasing level of aridity. The panels in the right column show 
the exact geographical location of the grid points used in the analysis. Wavelet cross¬ 
correlations T- in Figure 6 represent the coupling between T and 6 at a later time, so 
that T drives 6 if is signihcant at negative lags and vice versa, 6 drives T for positive 
lags. The t's^t is estimated based on a Paul kernel with m = 4, as an ensemble metric 
across the available range of years (1980-2011). The resulting correlation patterns can 
therefore be interpreted as ensemble averages of the 6-T coupling across scales. We use 
here the Paul wavelet based on its superior performance in resolving temporal variability. 
However, substantially similar results were obtained by using the Morlet kernel. At each 
location in space, rg^r is computed for the entire annual cycle (left panels), only Boreal 
summer (MJJA, central panels) and only Boreal winter (NDJF, right panels). Before 
inferring time-series of both T and 6 are normalized to have zero mean and unit 
variance. Therefore, harmonic and persistent oscillations are preserved, i.e. we are not 
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specifically analyzing seasonally de-trended anomalies, as one of the main goals is in fact 
to isolate the relative role of different harmonic components. 

As the expected correlation between T and 6 is mostly negative, the minimnm corre¬ 
lation {fmin) and its conhdence interval can represent a proxy of temporal asymmetry in 
the multiscale coupling also in this case. The conhdence interval of the ensemble seasonal 
Ts^t (winter or summer) are computed by considering each summer (or winter) as single 
time series, resulting in 24 annual realizations. For the full time series in contrast, the 
conhdence interval is calculated by using a sliding window approach. At each step, 
is estimated over a time-window of 5 years, sliding forward with a 1 year time step. The 
maximum scale at which Ts^t can be inferred in the annual and seasonal cases is limited 
by the size of each sample. Therefore we do not consider any periodicity larger than three 
months for seasonal samples since this modes are only partially sampled in the data. 

Panels a-c show for a location in central Sahara for the full time series (a), MJJA 
(b), and NDJF (c). Here, the extreme dryness of the soil is expected to inhibit coupling; 
consequently, a signihcant correlation can only be found for T driving 9 at the annual 
scale (circa 360 days) (see Figure 6a). This coupling, however, is mostly related to the 
seasonal cycle of temperature since in a hot desert climate, soil moisture retrievals mostly 
resemble a white noise signal and sensible heat dominates the exchange of energy between 
land and atmosphere. The extremely erratic behavior of the soil water content is captured 
by the symmetry of the minimum (absence of directionality in the coupling) and by 
the large variability in its conhdence interval (Figure 6a). Therefore, hot and dry desert 
climates can be seen here as a test for null directional coupling across scales, similarly to 
the synthetic case study reported in Figure 3c-d. 
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Correlation patterns radically change when we move to a location within the Sahelian 
sub-region of Mali (Figure 6d-f). In this case most of the precipitation (200-400 mm) falls 
during the summer months (MJJA), but it is sufficient to trigger the coupling between 
soil moisture and temperature across a wide range of temporal scales spanning a few days 
to months (Figure 6d,e). In contrast, the coupling disappears during the winter months 
(Figure 6f) due to soil desiccation. The seasonally dry character of this climate is captured 
by the second harmonic at half a year (around 180 days). During the entire year (6d), T 
leads 6 at time scales of one to 6 months, and 6 leads T at longer time scales. An analogous 
pattern can also be identihed in Northwestern Australia (Figure 6g-i), another seasonally 
dry location with a similar separation between humid and water-limited regimes. During 
Austral summer (NDJF), signihcant correlations are still conhned to the negative lags 
half-quadrant - i.e. in average, T drives 6 through evaporation - although the small scales 
(few days to 1 month) are clearly the most variable in terms of asymmetry of the coupling. 

This fact could be connected with a higher occurrence of 6 feedbacks on T at these 
scales throughout the twenty-one years of the ensemble (Figure 6e). Compared to Mali 
(Figure 6d), the variability of fmm, and consequently the uncertainty on directionality, 
is here less pronounced, i.e. the members in the ensemble (g and i) follow very similar 
correlation patterns at both small and large scales. During the rainy season (Austral 
winter. Figure 6h) 6 is not limiting evaporation, while during the summer the coupling 
extends throughout a wide range of scales (Figure 6i), with signihcant correlations also at 
positive lags {9 leading T) and at the hnest scales [Dirmeyer, 2011; Miralles et ai, 2012]. 

The separation between water-limited and wet regime is also present in the wavelet 
cross-correlation patterns of a temperate mid-latitude location in central France (Figure 
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6j-l). During all year (6j), the coupling is present at larger scales (> 90 days), and stronger 
for T leading 6 than 6 leading T. During the summer months (6k), negative correlations 
are found for negative lags, similar to the results for Mali (6i and 6e) and indicating 
an influence of T on 0 at all timescales. On the other hand, during the winter (61), 6 
leads T at large scales. The conhdence interval or fmin, however, mostly falls into null 
correlation regions and thus points to some uncertainties in these results. In general, the 
seasonal separation is weaker than in Mali (6e and f) and Northwestern Australia (6h and 
1). Finally, a site located not far from the coast of Gabon is displayed in Fignre 6m-o. The 
site is characterized by a tropical climate and strong moistnre advection from the Ocean, 
hence 6 does not represent a limiting factor in this case, and the only (positive) observed 
correlation is the one derived from the synchronization between the seasonal cycle of 
temperature and precipitation in the region [see e.g., Zhou et ai, 2014]. Overall the wavelet 
cross-correlation is able to captnre the scale-by-scale strength and directionality of the 
9-T conpling across different climatic regimes in a consistent way. Moreover, it allows the 
separation of the local scale contribution from the seasonal signal - often dominant in time 
domain statistics. Correlation results are stronger in transitional regimes, consistent with 
previous studies [Koster et ai, 2004; Seneviratne and Roster^ 2012; Miralles et ai, 2012], 
and interestingly, the d-T coupling propagates across all the analyzed scales dnring the 
warm season, when land-atmospheric interactions may be critical for high temperature 
extremes [Seneviratne et ai, 2006; Quesada et ai, 2012; Miralles et ai, 2014]. 

5. Conclusions 

We introduced a novel methodology to infer mnltiscale interactions from observations 
of dynamical systems that evolve over diverse temporal scales. The here adopted metric 
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- the wavelet cross-correlation - is based on the direct estimation of scale-by-scale 
correlations in the wavelet domain. The ability of to infer interactions across scales 

- and their directionality - was tested on different synthetic coupled systems and on a 
real-world case study of land-atmosphere interaction the coupling/feedback between soil 
moisture and near-surface air temperature. When applied to bi-variate auto-regressive 
vector models of increasing complexity the shows to be able to correctly reproduce 
the underlying directionality of the coupling at different temporal scales, and to distinguish 
fast/slow dynamic components within the simulated systems. In this context the term 
directionality is mainly used to indicate some sort of temporally lagged coupling at the 
considered scale, without any assumption on the causal structure of the observed process. 

However, it is clear that the ability of decomposing a coupling in its scale-by-scale com¬ 
ponents is an attractive feature of wavelet cross-correlation and can be used in disentangle 
the role of different dynamical processes in coupled geophysical systems. Besides direc¬ 
tionality is here mainly a synonymous for temporal asymmetry, and connections between 
causality and predictability ~ like in the case of more proper causality metrics like Granger 
causality [Granger, 1969] - were not yet explored. 

The application of wavelet cross-correlations to soil moisture and near surface air tem¬ 
perature shows interesting insights into the interaction of these two variables at different 
climate regimes. An interesting feature of this interaction, when observed through the 
lens of a multiscale correlation metrics like rg^r, is that the coupling between soil mois¬ 
ture and temperature, in the passage to a water limited regime, is active through an 
extremely wide range of scales, and that this feature represents a common signature of 
the coupling (and possible feedback) across very diverse climatic regimes. An interesting 
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question arising from this evidence could be whether this activation of the coupling across 
scales is a feature reproducible in climate models. A word of caution is hnally in order 
when we consider couplings taking place at large temporal scales (> 1 year), since they 
could be partially affected by the redundancy - and consequent auto-correlation effects - 
of the continuous wavelet transform. Overall results of this study demonstrate the poten¬ 
tial of wavelet cross-correlations to unravel the relationships between two environmental 
and climatic variables from a purely statistical perspective. The method here described 
can, in principle, be applied to observations from any region of the world, and to study 
soil moisture-temperature coupling or any other multivariate interaction across multiple 
scales. 
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Figure 1. Conceptual representation of wavelet coefficients Wi{t, sq) for i = X,Y, as & function 


of time t and a fixed sample scale sq (a-c), and wavelet cross-correlation (d-f) for X driving Y 


(a,d), instantaneous coupling (b,e) and Y forcing X (c,f) for different scales s and time-lags r 
(modified from Molini et al. [2010a]). The smallest temporal scale extracted in (d-f) is Si = ^, 
corresponding to the Nyquist frequency fu [see Torrence and Compo, 1998; Mallat, 2008; Molini 
et al, 2010a, for further details]. Forcing direction is taken homogeneous across scales. 
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Figure 2. Sample realization of the VAR(l) 
Lorentzian spectra P{f) (b) for Ci = —0.4 and 



system in equation (20) (a), and corresponding 
^2 = 0 
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Figure 3. Ensemble and “single-realization” estimates of Ts^t for the VAR(l) model in equa¬ 
tion (20) and different coupling parameters Ci and € 2 - Panels a and b respectively depict the 
ensemble estimates of Ts^t for a; —)■ j/ {Ci = 0, (72 = —0.9) and y ^ x {Ci = —0.4, (72 = 0). 
Regions with below the a = 99% signihcance level are masked in white. The same panels 
also show fmin (black empty circles) as defined in section 3.1 and the corresponding 99% conh- 
dence intervals (blue and red solid lines) as a function of scale. Wavelet correlation patterns for 
the uncoupled system {Ci = C 2 = 0) inferred from both a single realization, and an ensemble 
of simulations are shown in panels c and d, whereas panel e illustrates the sensitivity of f^m to 
the coupling strength and directionality for sample scales Si = 10 (red dots) and S2 = 120 (blue 
diamond) samples. Bottom and top x-axes in (e) respectively represent the coupling strength 
from from x to y {C 2 ) and y to x ((7i), whilst the red and blue bars are the conhdence intervals 
of fmin for s = Si and s = S 2 as in (a-b) and (d). Note how conhdence intervals tend to widen 
with the weakening of the coupling, and moving from hne to large scales consistently with panels 


a-b and d. 
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Figure 4. Ensemble estimates of Ts^t and fmm (same as in (a-b,d) of Figure 3) from the 
VAR(2) model in equation (21), obtained after pass-band filtering with the Morlet (a) and Paul 
(b) mother wavelets. Panel c shows the corresponding theoretical power spectra P{f) of the 
auto-regressive sub-spaces of x and y. Here, both the Morlet and the Paul wavelets are able to 
capture the presence of the pronounced periodicity at 1/5 of of the cycle, in good agreement with 
the coupling peak shown in panel c for the theoretical spectra. However, Morlet’s wavelet, being 
more localized, can better resolve the coupling peak and its directionality, at the expense of a 
higher redundancy at large temporal scales. 
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Figure 5. Same as Figure 4, but for a system with coupling and feedback occurring at 1/15 of 
the cycle and 1/6 of the cycle respectively, as described in equation (22). 
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Figure 6. Ensemble estimates of and f^m between soil moistnre 6 and air temperatnre 
T for five different geographical locations and climatic regimes. Right colnmn shows the exact 
geographical location of the grid points nsed in the analysis. Left colnmn reports the ensemble 
for the fnll time series at the specified location, while the central colnmns refers to Boreal 
snmmer (MJJA) and winter (NDJF) respectively. From top to bottom, panels are ordered by 
decreasing aridity. 
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